Two-fluid model for a rotating trapped Fermi gas in the BCS phase 
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We investigate the dynamical properties of a superfluid gas of trapped fermionic atoms in the 
BCS phase. As a simple example we consider the reaction of the gas to a slow rotation of the trap. 
It is shown that the currents generated by the rotation can be understood within a two-fluid model 
similar to the one used in the theory of superconductors, but with a position dependent ratio of 
normal and superfluid densities. The rather general result of this paper is that already at very low 
temperatures, far below the critical one, an important normal-fluid component appears in the outer 
regions of the gas. This renders the experimental observation of superfluidity effects more difficult 
and indicates that reliable theoretical predictions concerning other dynamical properties, like the 
frequencies of collective modes, can only be made by taking into account temperature effects. 

PACS numbers: 03.75.Kk,03.75.Ss,67.40.Bz 
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In the last few months, experiments with trapped 
fermionic 6 Li atoms made great progress. The fact that 
by using the Feshbach resonance the Fermi gas can be 
transformed into a Bose-Einstein condensate (BEC) of 
molecules, which can be cooled by evaporative cooling 
and afterwards transformed back into a Fermi gas, al- 
lows to reach extremely low temperatures of the order 
of 0.03T F 1], where T F = k B e F is the Fermi tempera- 
ture. This allows, among other things, a detailed study 
of the BEC-BCS crossover. In particular, a temperature 
of 0.03Tp should even be low enough to realize the BCS 
phase which is characterized by the condition A <C £.f, 
where A denotes the pairing gap. 

However, at present it is not very clear how the transi- 
tion to the BCS phase could be detected. While several 
observables related to collective oscillations (e.g. breath- 
ing modes) of the system have been investigated Q, , 
the most unambigious signatures of the superfluid BCS 
phase seem to be those which concern the rotational 
properties of the system Q. For instance, the moment of 
inertia of a slowly rotating Fermi gas was proposed to be 
a suitable observable for the detection of the BCS tran- 
sition p|. It should be mentioned that at present most 
of the theoretical predictions concerning possible experi- 
mental signatures of the BCS phase (e.g., 0, IE IE 00) 
neglect temperature effects as well as possible deviations 
from hydrodynamic behaviour due to the discrete level 
spectrum in the trap. 

In a previous article |E| we calculated the moment of 
inertia of a superfluid atomic Fermi gas in a slowly rotat- 
ing trap at finite temperature. There it turned out that 
the irrotational flow, which is characteristic for superflu- 
idity, is realized only in the limit when the gap A is very 
large compared with the temperature T and the level 
spacing hui of the trap. In all other cases, the velocity 
field has both, rotational and irrotational components. 
For example, if the level spacing hio is comparable with 
A, the current has a strong rotational component even 
at zero temperature. On the other hand, at nonvanishing 
temperature T, a certain fraction of the Cooper pairs is 
broken by thermal excitations. This leads to the well- 
known effect that the system behaves like a mixture of 



normal and superfluid components [Tfl ITU IT^ . Under 
rotation, the former behaves like a rigid body, while the 
latter can only have an irrotational velocity field. 

However, in the calculation of Ref. the gap A(r) has 
been replaced by a constant A corresponding to the aver- 
age diagonal matrix element of A(r) at the Fermi surface. 
While this averaging procedure seems to be justified in 
cases where only one oscillator shell participates in the 
pairing (intrashell pairing, A < fku), it is not well suited 
for the strong pairing regime (A > Hlj), where the prop- 
erties of the system can be described locally and depend 
on r via the spatial dependence of A(r) [13j. In partic- 
ular, the normal and superfluid fractions of the density, 
Pn/p and p s /p, should depend on r. To our knowledge 
this fact has not been taken into account in the existing 
literature. 

In this article, we will concentrate on the h — -> limit, 
i.e., we will neglect the quantum effect which is respons- 
able for the rotational component of the velocity field at 
zero temperature. Anyway, if the system is sufficiently 
large and if the temperature is not extremely low, this 
quantum effect becomes much smaller than the effect re- 
sulting from the thermally created normal component of 
the system. The important point is that we will now take 
into account the r dependence of the gap. In addition, 
we will not rely on the simplification made in our pre- 
vious work that the full potential (trap + mean field) is 
approximately harmonic. 

Let us briefly summarize the most important formulas 
(for more explanations and details see Ref. 0). We as- 
sume that equal numbers of atoms in two spin states are 
trapped in a harmonic potential, 



trap 



(i) 



The cigar-shaped form of the traps used in current ex- 
periments corresponds to lu z < = u> y . However, in 
order to force the system to rotate around the long axis, 
one has to break the axial symmetry, e.g., by using a 
rotating laser beam as "spoon" . We will model this by 
taking lu x ^ lu v . The mean- field single-particle hamilto- 
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nian minus the chemical potential reads 



h ° = 2^ + Vtra P^ + 9P( r ) - P ■ 



(2) 



the coupling constant g — A~Kh 2 a/m being proportional 
to the atom-atom scattering length a < 0, and p(r) be- 
ing the density per spin state. The order parameter in 
equilibrium is denoted Ao(r). 

Now we want to describe what happens if the trap 
is slowly rotating around the z axis with a rotation fre- 
quency f2 = f2e z . This is most easily done in the rotating 
reference frame. Then we still have a static problem, but 
the hamiltonian receives the additional term 



hi = —ttL z = — (fl x r) • p , 



(3) 



which we will treat as a small perturbation. Because 
of this term, the order parameter A(r) receives a phase 
exp[— 2i0(r)]. The explicit form of </>(r) is unknown for 
the moment and will be determined below. It is con- 
venient to eliminate this phase by a gauge transforma- 
tion, multiplying all single-particle wave functions by 
exp[i(j)(r)]. In this way the gauge transformed gap A(r) 
stays real, which in the case of a slow rotation implies 
that A does not change at all |U, i.e., A = A . On 
the other hand, this gauge transformation changes the 
momentum operator according to 



p = p h"V(f)(r) . 



(4) 



Hence, the price to pay for the real gap is an additional 
term in the hamiltonian. To linear order in the rotation 
frequency the new perturbation hamiltonian reads 



hi 



-VL Z - — (p • [V<Hr)] + [V0(r)] ■ p) • (5) 



In order to describe the system semiclassically, we 
make use of the Wigner-Kirkwood expansion. To that 

end we denote the Wigner transforms of ho, hi, etc., by 
h (r,p), hi(r, p), etc. We need also the Wigner trans- 
forms of the normal and abnormal density matrices in 
equilibrium, p (r,p) and Ko(r, p), as well as their devi- 
ations from equilibrium, pi(r, p) and Ki(r,p). (For the 
sake of brevity, we will occasionally omit the arguments r 
and p if there is no risk of confusion.) Furthermore we in- 
troduce the Poisson bracket of two phase-space functions 

z— x,y ,z 

Using the notations defined above, the terms linear in 
of the Hartree-Fock-Bogoliubov (HFB) or Bogoliubov-de 
Gennes equations up to linear order in H can be written 
as 

ih{h ,pi} + 2A «i = -ih{hi,p ) (7) 
ih{A 0l pi} - 2h ki — ih{hi, k } . (8) 



These are exactly the Eqs. (84) and (85) of Ref. [g. The 
main point of the present article concerns the solution of 
this system of equations in the case of an r-dependent 
gap A (r). 

First we eliminate hi by multiplying Eq. 10 by ho and 
Eq. ijHJ) by A and adding up the two resulting equations. 
Using the chain and product rules of differentiation, we 
then obtain 
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{£ 2 ,pi} = -h {hi,p } + A {hi,K } 



(9) 



where E 2 (r, p) = /ig(r,p) + Ag(r). To proceed further, 
we express po and kq in terms of ho and Ao- Within 
the Thomas-Fermi (TF) or local-density approximation, 
these relations read 

Po(r, P) = \- (1 - 2f[E(v, p)]) , (10) 



«o(r,p) 



Ao(r) 
2£(r,p) 



(l-2/[£(r,p)]), (11) 



where f(E) = 1 / '[exp(E / k B T) + 1] denotes the Fermi 
function. Although Eqs. i|l(J|) and (|ll|l are the solutions 
of the h — ► limit of the HFB equations, they are valid up 
to linear order in h 01 an d therefore consistent with Eqs. 
and ©. Inserting Eqs. (JTOJ and Jill) into Eq. we 
obtain, again after repeated use of chain and product 
rules of differentiation, the following simple equation: 

j{*.A>-5{*.*^}- < 12 > 

It is evident that this equation is solved by 

h = h%. (13) 

However, before this solution can be used, the gauge 
transformation which has been introduced in order to 
make the gap real must be inverted: 



P(r,p) = p[r, P + W0(r)] . 
To linear order in fi, this can also be written as 

P(r, p) =po[r, P + hV(f>(r)} + p x (r, p) 
=Po[r,P + hV<f>{r)} 



-(Qxr+-V^(r)).pi 

m ) dE 



E(r,p) 



(14) 



(15) 



(16) 



The last line has been obtained with the help of Eqs. (0 
and JT3J). 

The next step consists in calculating the corresponding 
current density per spin state, 



j(r) 



d 3 p p 

— p{ r ,P) 



(2nh) 3 m 



(17) 
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Using the explicit expression for p(r, p) given above, one 
easily obtains 



j(r) =-po(r) — V^(r) 



ft x r + - V0(r) ) / 

m /J 



dp 4 df 



6n 2 h 3P dE 



with 



PoO) 



d 3 p 
(2irhy 



Po(r,p) • 



(18) 



E(r,p) 

(19) 



(20) 



The current density can therefore be written in a more 
suggestive way in terms of normal and superfluid densi- 
ties p n (r) and p s (r), 



j(r) = p n (r)n x r - p s (r) — V<p(r) , 

m 



(21) 



if the normal and superfluid densities are defined accord- 
ing to the textbook result (Ref. |l5|. p. 459) as 



p n (r) = pq(t)-p s (t) = - 



dp 4 df 



6n 2 h 3P dE 



E(r,p) 



(22) 



In the BCS limit, i.e., if Ao(r) <C cf(t), where £ir(r) = 
P — V* rop (r) — 5Po( r ) denotes the local Fermi energy, the 
ratio p n (r) / po(r) becomes a function of only one dimen- 
sionless argument, T/T c (r), where T c (r) = 0.57A (r; T = 
0) denotes the local critical temperature (the existence of 
a local critical temperature is an artifact of the TF ap- 
proximation). This function, as well as the temperature 
dependence of the ratio A (r)/A (r;T = 0), is shown in 
Fig.ffl 

Up to now, the phase 0(r), which determines the veloc- 
ity of the superfluid component, is completely unknown. 
In Ref. |9j, this phase was determined by calculating Ai 
from ki and imposing the condition Ai = 0. Then it 
was shown that with this choice the continuity equation 
for the current was satisfied. Here we will adopt another 
method which is commonly used in the literature [Tlj and 
which consists in using the continuity equation. In the 
rotating frame the latter reads 



V-j(r) + p(r)-(nxr)-Vp(r) = 0, 



(23) 



where p(r) = in our case of a stationary rotation and 
p(r) = po(r) up to linear order in fi. Taking the di- 
vergence of Eq. 121fl . one can see that the normal fluid 
component drops out and it remains a continuity equa- 
tion for the superfluid component: 

-— V • [p,(r)V0(r)] - (n X r) • Vp s (r) = . (24) 

771 

In the case of a deformed harmonic trapping poten- 
tial, this equation can be solved analytically. To see 
this, remember that within the TF approximation the 
density po(v) and the gap Ao(r), and consequently also 
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FIG. 1: (a) Temperature dependence of the gap normalized 
to its value at zero temperature (A(0) = 1.76T C ). (b) Tem- 
perature dependence of the normal (solid line) and superfluid 
(deshed line) fractions of the superfluid system in the limit 
of weak pairing (A <C ej?)- These relations hold locally, if 
T c is interpreted as the local critical temperature, defined by 
T c (r) =0.57A(r;T = 0). 



the superfluid density p s , depend on r only via the local 
chemical potential pi oc (r) = p — Vtnip(r), i.e., p s (r) = 
p s [p — Vf rap (r)]. Hence Eq. (|24|) can be written as 



dp s 



dpio 



\VV trap (v)}- -V0(r) 
- -p,(r)V 2 
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ft x r 



(r) = 0. (25) 



In the special case of the harmonic potential , it can 
readily be shown that this equation has the same solution 
as in the simple case of constant Aq studied in Ref. [9j, 



(r) 
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(26) 



since this solution is independent of the form of p s (p)- 
The current density per spin state, j(r), is therefore given 
by 



j(r) = p n (r)n x r - p s (r)- 



• nV(r x r y ) . (27) 



Due to fact that po( r ) and Ao(r) [and consequently 
p s (r) and p n (r)] depend on r only via the local chemi- 
cal potential, it is sufficient to perform the TF calcula- 
tion for a spherical trap with the geometrically averaged 
trapping frequency Q — (uJxUjyUJz) 1 ^ 3 . In this spherical 
trap, of course, the density poi g a P Aq, etc., depend only 
on the distance from the center, i.e., po(r) = po(r), etc. 
(quantities related to the spherical trap will be marked 
by a bar). The corresponding quantities in the deformed 
trap can then be obtained from 



Po(r) = po 



— (— a> 2 r 2 + cu z r- 

| y I ^ y l 



(28) 
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FIG. 2: (a) gap Ao(r), (b) superfluid density p s (r), and (c) 
normal- fluid density p n (r) as a function of the distance r from 
a spherical trap with frequency Q for several temperatures: 
T = 0, 2, 4, 6 to/fc_B- For convenience, all quantities are given 
in harmonic oscillator units, i.e., r in Z^o = \J h/muo, Ao in 
to, p s and p„ in Zrf. The number of atoms is 36000, their 
ineraction strength is set to g = —fPlho/m. This choice of 
parameters allows us to compare the TF results (dashed lines) 
with results obtained from a numerical solution of the full 
HFB equations (solid lines). Note, however, that within the 
TF approximation one obtains qualitatively similar results in 
the case of more realistic parameters, if the temperatures are 
scaled with the corresponding critical temperature. 



and analogously for Ao(r), p s ( r )j etc. Note, however, 
that all this is true only within the TF approximation. 

In Fig. we show the normal and superfluid densities 
per spin state [p n and p s ) and the gap Ao in the spher- 
ical trap for different temperatures. The dashed lines 
correspond to the TF (h — > 0) results, while for the solid 
lines the gap Ao has been obtained by solving numeri- 
cally the HFB equations ^(|. At zero temperature, the 
trapped Fermi gas is completely superfluid, i.e., p s = pQ 
and p n = 0. At low but non-vanishing temperature 
(T = 2huj ps 0.36T C , T c w 5.5to being defined as the crit- 
ical temperature within HFB), a normal fluid component 
appears near the surface, since there the gap is smallest 
and consequently the Cooper pairs are most easily bro- 
ken by thermal excitations. If the temperature increases 
further (T = 4fcl> k 0.72T C ), the normal-fluid compo- 
nent starts to extend over the whole volume, and finally, 
slightly above the critical temperature (T = 6Huj) the 
superfluid component vanishes completely (solid lines). 
The small superfluid region near the center of the trap 
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FIG. 3: Current density divided by the angular velocity of 
the rotation, j(r)/fl, in the xy plane (r z = 0) of a trap with 
uj x /ujy — 0.7 and LUz/uJr = 0.03 (uj r — y/ui x ui y ) at four different 
temperatures. The remaining parameters are the same as in 
Fig. HThe length of the arrow displayed in the lower right 
corner of each figure corresponds to a current density of j — 
20fl/lf lo . The coordinates r x and r y are given in units of Iho 
and the temperatures in units of to/fcs. 



which survives at this temperature within the TF ap- 
proximation (dashed line) is a consequence of the fact 
that the critical temperature predicted by the TF ap- 
proximation is higher than the HFB one |l3 - However, 
apart from this point, one can say that in general for p n 
and p s the agreement between TF and HFB is reason- 
able and better than for the gap Ao itself. The reason 
why the agreement between TF and HFB is better for p n 
and p s than for Ao is that near the critical temperature 
the temperature dependence of p a j p is much weaker than 
that of A/A(T = 0), cf. Fig.HJ 

In order to make the comparison between the gaps Ao 
calculated within the TF approximation and by solving 
the full HFB equations, we had to choose a rather mod- 
erate number of particles, N — 36000, for which the HFB 
calculation is feasible. If we had not been interested in 
the comparison between the results obtained with the 
HFB and TF gaps, we could of course have shown the 
TF results for arbitrarily large numbers of particles. We 
emphasize that even for much larger numbers of particles 
the qualitative behaviour of the TF results remains un- 
changed, provided that the coupling constant g is tuned 
such that the condition A <§C ep (BCS condition) re- 
mains satisfied and the temperatures are scaled with 
respect to the critical temperature. 

Using the spherical density profiles and Eqs. and 
(|27|l . we can immediately calculate the current distri- 
bution j (r) . In Fig. [3] we show the current in the xy 
plane (r z = 0) for a deviation from axial symmetry of 
= 0.7 at several temperatures. (For the cases 
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FIG. 4: Moment of inertia, O of a trapped Fermi gas as 
function of temperature. The parameters are the same as in 
Figs. [5] and [3] The moment of inertia and the temperature 
are given in harmonic oscillator units, i.e., in rrd\ and T 
in to/fc_B. The transition temperature to the BCS phase lies 
at approximately 5.5to. The solid line has been obtained by 
performing a full HFB calculation for po(r) and Ao(r), while 
for the dashed line the TF results have been used. 



shown, the current densities obtained from the HFB and 
TF density profiles are indistinguishable within the res- 
olution of the plot.) At zero temperature, the current 
is irrotational and rather weak (it vanishes in the limit 
of axial symmetry, u x = u> y ). In the surface region the 
current reaches its ordinary (rigid-body) form already at 
T = 1.5to » 0.27T C ). At T = ihui rj 0.54T c the current 
shows almost everywhere the rigid-body behavior, only 
near the center it is still a little bit weaker than in the 
normal phase, T — 6to. 

Let us now look at the temperature dependence of the 
moment of inertia, O, which is defined as 



9 



(4) 
ft 



2m 
"ft" 



d 3 r[r x j y (r) -r v j x (r)] 



(29) 



The factor of two is a consequence of our convention that 
j denotes the current density per spin state. Using again 
Eqs. I|28|l and l|27|) . we can express the moment of iner- 
tia in terms of the density profile in the corresponding 
spherical trap as a simple radial integral: 



87TTO/U/ 2 

~3~ w 



or 



dr r 



Pn(r) 



+ Gzri)^4 (30) 



\UJ A j, + OJ, 

In Fig. 0] we show the moment of inertia for the same 
set of parameters that were already used in Figs. |21 and 
OH as a function of temperature. The solid line has been 
calculated by using the HFB density profiles, while the 
dashed line was obtained from the density profile within 
TF approximation. One can see that the moment of in- 
ertia decreases strongly as the temperature goes to zero. 
The limiting value at zero temperature is determined by 
the deformation of the trap in the xy plane, 



6(T = 0) 



e 



rigid ■ 



(31) 



where & r igid denotes the corresponding rigid-body mo- 
ment of inertia [which can be obtained from Eq. (I3UI) by 
putting p n — po and p s = 0]. In our case of oj x /ll> v = 0.7, 
we have 0(T = 0) ~ 0.12O rtg id. An important point to 
notice is that, coming from high temperatures, one does 
not observe an appreciable change of the moment of iner- 
tia until one reaches temperatures far below the critical 
one. The reason for this effect is that the main contri- 
bution to the moment of inertia comes from the outer 
regions of the trapped gas, where the order parameter 
becomes small and where consequently the normal-fluid 
fraction is large even far below T c . The discrepancy be- 
tween the HFB (solid line) and TF results (dashed line) 
below w 2hw can be traced back to the effect that within 
the TF approximation the gap near the surface vanishes 
already at very low temperature, such that the normal- 
fluid fraction near the surface is overestimated within 
TF. 

To conclude, we have applied the two-fluid model 
known from the theory of superconductivity [ToL ITU Il2]| 
to the case of ultracold trapped fermionic atoms in the 
BCS phase. In contrast to the usual situation, the ra- 
tio of the normal and superfluid densities is explicitly 
position dependent due to the inhomogeneous trapping 
potential. Specializing to the case of a slowly rotating 
system, we have shown that the linear order in h of the 
linear response equations gives a current which can be 
decomposed in a natural way into normal and superfluid 
components. The normal component appears as a con- 
sequence of Cooper pairs which are broken by thermal 
excitations already below the critical temperature T c . 
We have shown that especially the outer region of the 
trapped gas behaves essentially as if it was normal-fluid, 
even far below T c . Only the central region of the gas 
keeps its superfluid character up to T c . 

As a consequence, the moment of inertia decreases 
more slowly as it was previously expected if the tem- 
perature is lowered below T c , i.e., the effects of superflu- 
idity become visible only far below T c . This important 
but in a certain sense negative result will apply analo- 
gously to other observables which are mainly sensitive 
to the surface of the system, like, e.g., collective modes. 
For example, the theory presented here was used in Ref. 
|18| in order to explain the temperature dependence of 
the strength of the response function for the so-called 
"twist mode" . There the effect was even more dramatic, 
since the relevant integral contained an r 6 weight factor 
instead of r 4 in Eq. I|30fl . 

We are therefore convinced that it is not justified to 
compare the experimentally measured frequencies of col- 
lective modes directly with theoretical predictions ob- 
tained for zero temperature, as it is done in the current 
literature [l|,l2l- We ra f ner expect that the temperature 
dependence is important and can be predicted by gener- 
alizing the two-fluid model presented here to the dynamic 
case, i.e., by performing the Wigner-Kirkwood expansion 
of the time-dependent HFB (TDHFB) equations up to 
linear order in H 0, . This leads to a generalization 
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of the Vlasov equation for the normal phase, which re- 
sults from the Wigner-Kirkwood expansion of the time- 
dependent Hartree-Fock (TDHF) equation up to linear 
order in h. 

However, one should keep in mind that the Thomas- 
Fermi approximation for the ground state as well as 
the Wigner-Kirkwood expansion of the dynamical equa- 
tions (i.e., the generalized Vlasov equation and the su- 
perfluid hydrodynamics to which it reduces in the zero- 
temperature limit) depend on the assumption fkUi <C A 
for i = x,y, z. Concerning the validity of the Thomas- 
Fermi approximation for the ground state, the condition 
huji <C A has been inferred from the requirement that 
the coherence length £ = hyp/itlS. (vf being the Fermi 
velocity) must be much smaller than the typical length 
scale of the system, which is approximately given by the 
Thomas- Fermi radius Rtf — y/2/i/mujf [13. However, 
it is less evident where the assumption Hu>i <C A enters 
into the description of the dynamics of the system within 
the generalized Vlasov equation. 

To give a specific example, in Ref. 0, quantum cor- 
rections to the moment of inertia of higher orders in 
H(oj x ± oj y )/A were discussed. Also in the case of the 
strength of the twist mode mentioned above, the fully 
quantum- mechanical ( "microscopic" ) calculation showed 
deviations from the two-fluid model, especially at very 
low temperatures. In both cases, the corrections act as if 
the normal- fluid component of the system was larger than 
predicted by Eq. i|22|) and in particular non-vanishing 
even at zero temperature. In a certain sense the accel- 



erations acting on the Cooper pairs during their motion 
through the inhomogeneous potential seem to have a sim- 
ilar pair-breaking effect as the thermal excitations which 
are responsable for the normal-fluid component given by 
Eq. 122(1 . From a completely different point of view, Eq. 
(|22|l is usually derived by looking at the long-wavelength 
limit (g£ <C 1) of the current-current correlation func- 
tion in a homogeneous system [12J, Il5| . In the trapped 
system, however, the wave vectors must be of the order 
q ~ 1/Rtf, and we recover the condition £ <C Rtf- 

Deviations from superfluid hydrodynamics (T — 0) or 
from the two- fluid model (T > 0), respectively, may 
therefore be especially important in the case of the 
strongly elongated traps used in current experiments, 
which have rather high radial trapping frequencies uj x 
and ujy . Therefore this kind of quantum effects should be 
studied in more detail. In the case of collective modes, 
this could be done, e.g., by comparing systematically 
the results obtained in quantum mechanical quasiparti- 
cle RPA (QRPA) calculations [2(j with those of hydro- 
dynamics. 
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